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I.  Introduction 


Over  the  last  two  decades,  the  Direct  Simulation  Monte  Carlo  (DSMC)  method  has  gradually  evolved  into 
a  powerful  and  versatile  technique  that  has  been  successfully  applied  to  many  complex  problems  of  modern 
rarefied  gas  dynamics,  primarily  related  to  high-altitude  aerodynamics  and  microscale  flows.1  However, 
nearly  all  current  applications  of  the  DSMC  method  deal  with  gas  flows  interacting  with  stationary  or 
pseudo-stationary  gas-surface  interfaces.  Two  exceptions  are  the  work  of  E.  Oran  et  al.,2  who  modeled  a 
rotating  ellipse  using  a  nronotonic  Lagrangian  grid  technique,  and  the  work  of  Gallis  et  al.,3  who  examined 
flow  over  moving  microbeams  using  both  continuum  (Navier-Stokes)  and  kinetic  (DSMC)  approaches.  There 
is  therefore  a  large  gap  between  the  relatively  mature  state  of  continuum  modeling  of  moving  solid  interfaces 
(such  as  fixed-wing  aircraft,  rotor  aircraft,  and  turbonrachinery  applications)  and  particle-based  approach 
such  as  DSMC,  where  capabilities  such  as  six  degrees  of  freedom  (6-DOF)  body  motions  and  conjugate  heat 
transfer  are  in  relatively  early  stages.  There  are  several  major  challenges  that  limit  current  application  of 
the  DSMC  method  to  study  gas  flows  with  moving  solid  boundaries.  These  include  algorithmic,  mostly 
grid-related  difficulties  of  efficient  realization  of  the  DSMC  method,  especially  for  a  free  type  body  motion 
with  arbitrary  shape  changes  and  numerical  problems  that  arise  due  to  the  subsonic  nature  of  these  flows. 

To  help  extend  the  DSMC  method  as  a  practical  tool  for  rarefied  time-dependent  flows,  we  present 
results  for  a  generic  high  altitude  (~100  km)  sounding  rocket  stage  separation  scenario,  where  the  relative 
atmosphere-rocket  speed  is  5  knr/s.  We  examine  two  cases:  a  liquid  propellant  with  25  kN  thrust  and  a  solid 
propellant  with  34  kN  thrust.  In  our  calculations,  the  upper  stage  thruster  fires  on  a  lower  stage  (0.75m 
diameter  3  m  long  cylinder)  initially  at  2  meters  separation  with  an  initial  relative  velocity  of  2  nreters/second. 
These  parameters  along  with  a  high  thruster-related  acceleration  represent  conditions  typical  for  a  very  fast 
stage  separation,  and  are  used  to  emphasize  the  transience  of  the  process.  Figure  1  shows  the  basic  stage 
separation  scenario. 

Among  the  primary  motivations  of  this  study  is  to  determine  whether  a  DSMC  approach  is  a  practical 
numerical  scheme  for  a  typical  stage  separation  in  a  rarefied  atmosphere,  and  to  test  how  well  quasi-steady 
state  treatments  compare  to  fully  time-dependent  results.  The  initial  stage  separation  velocities  and  other 
parameters  in  the  flow  are  therefore  chosen  to  stress  steady-state  assumptions,  but  at  the  same  time,  the 
parameters  reflect  practical  scenarios.  Another  important  element  in  the  present  study  is  the  determination 
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of  forces  and  heat  loads  on  the  lower  stage,  as  well  as  determination  of  the  radiation  signature  due  to 
interaction  of  the  upper  stage  plume  with  the  atmosphere  and  impingement  of  the  upper  stage  plume  on  the 
lower  stage.  We  have  therefore  performed  a  matrix  of  calculations  with  the  above  generic  stage  separation 
scenario  which  includes:  1)  steady  state  with  lower  stage,  2)  unsteady  with  moving  lower  stage,  3)  unsteady 
with  a  fixed  lower  stage,  and  4)  unsteady  with  no  lower  stage. 

One  of  the  chief  results  we  have  found  is  that  numerical  convergence  can  be  practically  achieved  within 
our  DSMC  approach  for  our  unsteady  and  steady  calculations.  Another  significant  finding  is  that  the  steady- 
state  runs  with  the  lower  stage  at  a  fixed  distance  (quasi-steady)  closely  match  the  full  unsteady  results. 
Because  the  scenario  conditions  are  at  the  upper  end  of  practical  applications  in  terms  of  velocities  and 
thrusts  and  would  be  most  likely  to  show  unsteady  effects  if  they  were  present,  our  results  show  that  quasi 
steady-state  calculations  are  adequate  to  describe  most  rarefied  stage-separation  scenarios.  This  makes 
examination  of  such  phenomenon  much  more  practical  with  DSMC.  The  major  features  of  the  flow  are 
a  plume-atmosphere  shock/mixing  layer,  the  plume-gas  impingement  shock  on  the  lower  stage,  and  the 
particle  flow  and  impingement  on  the  lower  stage.  The  plume-atmosphere  shock,  which  reaches  10,000 
K  in  translational  temperature,  appears  to  be  completely  independent  of  the  presence  lower-stage  for  the 
conditions  examined.  We  also  generated  radiation  maps  of  our  flow-fields  for  a  sensor  viewing  the  flow 
from  a  large  distance.  We  note  that  the  stage  separation  problem  has  been  treated  before  in  the  context  of 
continuum  fluid  dynamics.4  To  our  knowledge,  however,  the  present  work  is  the  first  application  of  DSMC 
to  time-dependent  stage  separation  in  a  rarefied  atmosphere. 

II.  Stage  Separation  Forces 

Before  discussing  the  particulars  of  the  numerical  approach,  it  is  useful  to  define  the  forces  and  other 
dynamical  variables  which  play  a  role  in  the  stage  separation  scenario  being  considered.  Prior  to  stage 
separation,  the  rocket  is  oriented  in  the  es  direction.  The  lower  stage  turns  off,  causing  the  vehicle  to  enter 
free  fall  under  gravitational  acceleration  g.  At  this  point,  both  stages  are  attached  to  each  other.  Then  a 
separation  mechanism  provides  an  impulse  of  force,  causing  the  upper  and  lower  stages  to  separate  from 
each  other  with  initial  relative  velocity  vsep,  where  the  individual  stages’  velocities  are  V2lS  —  VijS  =  vsep. 
The  separated  stages  are  in  free  fall  for  time  duration  ts,  after  which  the  upper  stage  propulsion  fires. 
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Figure  1.  Free-body  diagram  of  the  two-stage  system  after  separation.  Subscript  2  denotes  the  upper  stage 
and  subscript  1  denotes  the  lower  stage. 


Figure  1  shows  the  free-body  diagram  for  the  system.  The  gas  flow  exerts  force  on  the  bodies  through 
thrust  (T),  plume  impingement  (Fp),  and  drag  (Fdi,  Fd2).  Gravity  of  course  exerts  a  force  proportional 
to  object  mass,  Fgi  =  m^g.  Then,  if  we  define  t  =  0  as  the  time  immediately  after  separation,  the  stage 
velocities  can  be  related  to  the  forces: 


vi  (t) 


v2(t) 


v1)S+  fFgl+F  P  +  Fdl  dt 

J  m  i 


f  Fp  Fdi 

vi,s  +  /  g  H - 1 - dt 

J  mi  mi 


V2,; 


Fg2  +  T  +  Fd2 
m2 


dt 


(1) 

(2) 

(3) 


t 


=  v2jS  + 


m2  m2 


dt 


(4) 


o 

We  assume  that  the  distance  between  stages  is  not  so  large  as  to  produce  a  discernible  change  in  the  gravity 
g  vector.  We  further  constrain  vi+2(t  =  0)  and  vsep  to  be  parallel  to  es,  so  that  there  is  no  angle-of-attack. 
(The  apparent  direction  of  the  drag  force  vector  in  Fig.  1  is  a  generalization.)  The  resulting  separation 


velocity  is  then  independent  of  gravity: 

a  r+\  r+\  i  /  T  +  Fd2  Fp  +  Fdi 

Av(t)  =  v2(i)  -  vi(t)  =  vsep  +  / - - -  dt  (5) 

J  m2  mi 

o 
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In  this  simple  linearly-translating  system,  we  may  use  a  gravitationally-accelerated  reference  frame.  Gravity’s 
main  effect  would  be  to  change  the  velocity  of  the  atmosphere  relative  to  the  vehicles,  but  for  short  time 
durations  (1  s)  the  error  is  less  then  10  m/s,  which  is  negligible  compared  to  the  original  velocity  of  5000  m/s. 

Because  we  are  at  100  km  altitude  and  the  time  under  consideration  is  relatively  short  (1  s),  the  atmo¬ 
spheric  drag  forces  can  be  neglected  compared  to  the  thrust  and  plume-induced  forces,  further  reducing  the 
time  dependent  separation  relation: 

t 

A v(i)  «  vsep  +  f  — - -  dt  (6) 

J  m2  m  i 

o 

Important  to  this  assumption  is  the  orientation  of  the  vehicle  axis  parallel  to  the  direction  of  travel.  Although 
the  drag  force  is  small  compared  to  the  thrust  and  plume-impingement  forces,  the  system  is  more  sensitive 
to  lateral  (i.e.,  _L  es)  perturbations  in  alignment  of  the  stages,  as  they  could  cause  the  lower  stage  to  tumble. 

III.  Numerical  Approach 

The  DSMC  method  is  used  to  model  the  stage  separation,  with  the  lower  stage  moving  away  from  the 
upper  stage,  and  a  number  of  changes  need  to  be  introduced  into  the  standard  DSMC  algorithm  to  account  for 
that  motion.  The  principal  changes  include  (i)  stage  motion  under  the  impact  of  the  two  main  contributing 
forces,  T  and  Fp,  (ii)  collisions  of  molecules  with  the  moving  stage,  and  (iii)  displacement  of  molecules  by 
the  moving  stage.  In  order  to  incorporate  these  changes,  the  following  algorithms  are  used. 

(i)  Stage  motion. 

The  stage  will  move  every  N  time  steps  (N  =  100  was  used  in  this  work),  with  the  condition  that  the 
stage  displacement  is  small  compared  to  the  distance  between  the  stages  and  lower  stage  size.  This  condition 
will  ensure  that  the  impact  of  discretization  of  the  stage  motion  on  the  flow  field  is  small.  The  distance  that 
the  stage  is  displaced  at  is  proportional  to  N At  and  the  relative  stage  velocity  Av  described  in  the  previous 
section.  Here,  At  is  the  time  step,  and  the  velocity  Av  depends  on  the  value  of  the  plume  force  Fp  and  the 
constant  thrust  force  T. 

Since  the  plume  force  Fp  changes  as  the  stage  separation  increases,  the  current  value  of  Fp  is  used  to 
estimate  the  plume  force  related  contribution  to  Av.  The  current  value  was  obtained  through  the  averaging 
of  the  surface  properties  using  gas-surface  collisions  that  occur  over  the  last  M  timesteps  before  the  actual 
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stage  motion  (M  =  10  was  used  in  this  work).  In  the  axisynnnetric  simulations  presented  below,  only  the 
drag  component  of  the  plume  force  is  accounted  for,  and  the  lower  stage  moves  only  along  the  plume  axis. 
In  the  two-dimensional  simulations,  both  lift  and  pitching  moment  are  included,  and  the  stage  is  allowed 
to  move  off  the  plume  axis  and  to  rotate.  This  kind  of  motion  could  occur  due  to  DSMC  statistical  errors, 
which  may  create  asymmetric  forcings.  Therefore  such  tests  are  a  measure  of  the  robustness  of  our  DSMC 
moving  body  algorithms.  For  the  rotation,  the  center  of  mass  was  assumed  at  the  geometric  center  of  the 
lower  stage.  In  the  algorithms  presented  below,  the  stage  rotation  is  not  accounted  for,  and  only  translational 
motion  is  assumed  to  be  contributing  to  gas-surface  collisions  and  molecule  displacement. 

(ii)  Gas-surface  collisions. 

Since  the  stage  motion  is  modeled  in  the  reference  frame  of  the  upper  stage,  gas-surface  collisions  with  the 
upper  stage  are  modeled  as  usual.  For  the  gas-surface  collisions  with  the  lower  stage,  molecular  velocities  after 
collision  (normal  and  tangential  components)  are  sampled  from  the  Maxwell  distribution  function  (diffuse 
reflection),  which  reflects  the  assumed  lower  stage  temperature,  plus  the  added  velocity  that  corresponds  to 
the  translational  velocity  of  the  body. 

(iii)  Molecule  displacement. 

To  displace  molecules  by  the  moving  stage,  at  the  end  of  each  N  timestep  the  following  simple  algorithm, 
which  conserves  mass  and  energy,  is  used. 

1.  For  each  molecule  in  the  computational  domain,  check  if  it  is  close  enough  to  the  lower  stage’s  center 
of  mass  (this  excludes  from  consideration  a  number  of  molecules  far  from  the  surface). 

2.  If  a  molecule  is  close  to  the  center  of  mass,  compute  its  collision  time  with  the  stage,  tc,  assuming  the 
stage  does  not  move  and  the  molecule  moves  with  the  velocity  (— Avx,0,0)  in  the  axisymmetric  case  and 
(— Avx,  —  Avy,  0)  in  the  2D  case.  Here,  (Avl;,  Av.y)  are  the  components  of  the  current  stage  velocity. 

3.  If  tc  <  NAt,  compute  new  coordinates  of  the  molecule  as  if  it  had  a  velocity  — Av  and  experienced  a 
specular  collision  with  the  surface,  as  —A vfc  +  Av(NAt  —  tc)  +  AvNAt  =  2Av(NAt  —  tc). 

To  implement  these  algorithms,  we  use  the  SMILE  DSMC  code,5  and  extend  it  with  the  free-cell  method,6 
which  avoids  a  potentially  numerically  intensive  and  complex  grid  reconstruction  step  during  the  stage 
motion. 


6  of  21 


American  Institute  of  Aeronautics  and  Astronautics 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


IV.  Flow  Conditions  and  Parameters  of  the  Approach 


We  consider  a  generic  sounding  rocket  that  consists  of  a  cone-cylinder  upper  stage  and  a  cylinder  lower 
stage.  The  diameter  of  both  stages  is  assumed  to  be  0.75  m,  the  length  of  the  upper  and  lower  stage  is  3  m, 
and  the  conical  part  is  1  m  long.  The  upper  and  lower  stage  masses  are  350  kg  and  125  kg,  respectively. 
This  results  in  a  thrust-based  acceleration  of  the  upper  stage  of  about  70  m/s 2 .  We  use  an  atmosphere  of 
N2  at  conditions  approximating  100  km  altitude  and  perform  calculations  up  to  1  second  when  the  lower 
stage  is  50  m  from  the  upper  stage.  The  staging  is  calculated  at  a  free  stream  velocity  of  5  km/s,  and  the 
atmosphere  is  assumed  to  have  a  temperature  of  190  K  and  number  density  of  8.5xl018  mol/m3.  The  stage 
surface  temperature  is  uniform  at  500  K.  Initial  stage  separation  distance  and  velocity  are  2  m  and  2  m/s, 
respectively. 

Two  types  of  thrusters  are  considered,  liquid  and  solid  propellant  thrusters.  In  both  thrusters,  the  gas 
stagnation  pressure  and  temperature  are  35  atm  and  3000  K,  respectively.  To  simplify  flow  analysis,  the 
gas  plume  was  presented  as  a  75%  CO  and  25%  H20  mixture.  The  nozzle  throat  and  exit  diameters  are 
7  cm  and  50  cm,  respectively,  and  the  diverging  part  half-angle  is  15  deg.  The  resulting  thrust  for  the  liquid 
propellant  thruster  is  about  25  kN.  In  the  solid  propellant  thruster,  in  addition  to  the  component  gas,  two 
alumina  particle  sizes  are  modeled  with  diameters  of  0.2  /jm  and  6  /zm.  The  aluminum  particle  loading  is 
about  20%,  and  the  resulting  thrust  is  about  34  kN. 

The  flow  inside  the  nozzle  is  modeled  using  CFD++  software.7  CFD++  is  a  flexible  computational 
fluid  dynamics  software  suite  for  the  solution  of  steady  and  unsteady,  compressible  and  incompressible 
Navier-Stokes  equations,  including  multi-species  capability  for  perfect  and  reacting  gases.  In  this  work,  a 
reacting-gas  compressible  Navier-Stokes  solver  is  used  with  second  order  spatial  discretization  and  implicit 
time  integration;  a  general  multi-phase  capability  has  been  used  that  provides  an  Eulerian  description  of  the 
disperse  phase  (particulates).  The  macroparameter  profiles  at  the  nozzle  exit  are  obtained  using  a  single¬ 
block  rectangular  grid  with  a  total  of  about  20,000  nodes.  The  inflow  and  outflow  boundary  conditions 
are  stagnation  pressure  and  temperature,  and  backpressure  imposition  (the  value  of  1  Pa  is  assumed), 
respectively.  A  symmetry  condition  is  used  at  the  plume  axis.  The  isothermal  boundary  with  a  surface 
temperature  of  1,000  K  is  assumed  at  the  nozzle  surface. 

The  plume  expanding  into  a  rarefied  atmosphere  of  100  km,  of  N2,  is  computed  using  SMILE  computa- 
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tional  tool5  based  on  the  DSMC  method.  A  starting  surface  at  the  nozzle  exit  plane  is  used  to  model  plume 
boundary  conditions  with  parameters  obtained  using  CFD++.  About  10  million  molecules  and  300,000 
background  collision  cells  were  used.  The  number  of  sampling  cells  is  about  350,000.  A  time  step  of  10-5 
seconds  was  used  and  the  body  movement  occurred  every  100  time  steps  (0.001  seconds).  The  change  in 
stage  separation  distance  over  this  time  is  smaller  than  0.1%,  which  is  considered  adequate  in  terms  of  nu¬ 
merical  accuracy.  Note  that  near  the  nozzle  exit  the  cell  size  and  the  time  step  do  not  satisfy  the  conventional 
DSMC  requirements,  and  a  collision  limiter8  was  used  there  in  order  to  reduce  computational  time.  To  check 
the  accuracy  of  results,  a  single  calculation  was  conducted  with  about  four  times  the  number  of  molecules 
as  above.  The  difference  between  the  results  was  found  to  be  within  statistical  scatter.  The  numerical 
approximations  used  in  our  DSMC  implementation  therefore  should  not  be  a  significant  source  of  error. 

The  VHS  model  is  used  for  gas-gas  collisions,  and  a  diffuse  reflection  model  is  used  for  gas  scattering 
with  surfaces.  The  gas-particulate  interaction  is  modeled  with  the  approach  of  Gallis  et  al.9  using  two-way 
coupling.10  The  DSMC  modeling  accounts  for  the  plume  impingement  forces  and  lower  stage  movement. 
Lacking  experimental  guidance  at  the  relative  velocities  of  interest,  we  assume  a  diffuse  reflectance  model 
for  particulate  reflection  from  the  lower  stage.  The  lack  of  a  validated  or  generally  accepted  particulate- 
surface  interaction  model  is  a  significant  gap  in  the  modeling  of  such  scenarios.  Most  calculations  used  axial 
symmetry.  All  axisymmetric  computations  have  been  conducted  using  radial  weights.  Species  weights  are 
used  to  reduce  the  statistical  scatter  for  the  free  stream  species  and  particulates. 

V.  Stage  Separation  Driven  by  a  Liquid  Propellant  Thruster 

The  first  set  of  computations  was  conducted  for  a  25  kN  liquid  propellant  thruster,  and  the  stage  sep¬ 
aration  was  modeled  over  the  first  second  after  thruster  ignition.  In  order  to  understand  the  role  of  lower 
stage  motion  on  the  flow  and  to  examine  several  key  modeling  assumptions,  the  following  four  cases  were 
considered.  First,  unsteady  full  stage  motion  is  modeled  as  a  result  of  the  thrust  force  T  propelling  the 
upper  stage  and  the  force  Fp  from  the  plume  impinging  on  the  lower  stage.  Second,  the  flow  was  computed 
for  an  unsteady  plume  development  but  a  fixed  body  geometry  (two  geometries  were  considered  for  the  time 
after  ignition  of  0.3  s  and  1  s).  Third,  an  unsteady  plume  development  was  simulated  at  0.3  s  and  1  s  with 
no  upper  and  lower  stage  geometries  involved  (plume  only  configuration).  Finally,  a  steady  state  flow  was 
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calculated  for  stage  separation  distances  of  7  m  and  42  m  that  correspond  to  the  above  two  time  moments  in 


the  moving  stage  configuration.  To  simplify  analysis,  the  lower  stage  was  moved  relative  to  the  upper  stage, 


and  the  coordinates  of  the  upper  stage  were  fixed  in  all  simulations. 
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Figure  2.  Liquid  propellant  thruster:  translational  temperature  (K)  at  0.3  s  after  thruster  ignition.  Here  and 
below,  the  direction  of  the  free  stream  velocity  is  from  left  to  right,  going  from  -x  to  +x. 


The  results  for  the  four  cases  are  given  in  Fig.  2,  where  the  translational  temperature  fields  calculated 
over  all  gas  species  are  presented  that  correspond  to  a  time  moment  of  0.3  s  after  ignition.  Only  the  upper 
half  of  the  symmetric  flow  is  shown,  and  the  plume  axis  and  both  stage  centerlines  are  aligned  along  the 
horizontal  X  axis.  The  leading  edge  of  the  upper  stage  is  at  X=0,  and  that  of  the  lower  stage  is  at  X«10  m, 
so  that  the  separation  between  them  is  about  7  m.  The  radius  of  the  two  stages  is  0.345  nr,  and  so  may  not  be 
visible  on  the  scale  of  the  figure.  The  temperature  flowfields  show  two  shock  regions:  the  plume-atmosphere 
mixing  layer  in  a  wide  swath  starting  at  the  upper  stage  nozzle  exit  and  the  plume-lower  stage  interaction 
region  starting  at  the  leading  edge  cap  of  the  lower  stage.  The  free  stream  velocity  (5  knr/s)  is  about  two 
times  larger  than  the  typical  plume  exhaust  velocity  (2.5  km/s),  which  helps  explain  the  higher  temperature 
in  the  first  region  compared  to  the  flow  near  and  after  the  lower  stage.  The  maximum  temperatures  in  these 
two  regions  are  about  7,200  K  and  2,300  K,  respectively.  The  plume-atmosphere  shock  mixing  layers  are 
very  similar  in  all  four  cases.  However,  the  mixing  layer  in  the  moving  stage  case  appears  to  be  thinner 
than  in  the  steady  state  case  primarily  due  to  the  small  sampling  size  and  not  the  actual  stage  separation 
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process.  Statistical  reasons  are  also  behind  a  small  difference  in  temperature  in  the  lower  stage  wake  flow, 
where  the  temperature  is  lowest  for  the  moving  stage,  and  highest  for  the  steady  state  case.  The  average 
number  of  particles  in  a  sampling  cell  in  the  wake  flow  region  was  on  the  order  or  less  than  ten,  which 
explains  the  inaccuracy  of  the  temperature  prediction.  The  five  times  longer  sampling  time  in  the  fixed 
body  case  results  in  higher  temperatures  in  the  wake  compared  to  the  moving  stage  (in  the  latter  case,  the 
temperature  calculation  was  meaningless  for  a  number  of  cells  where  less  than  two  molecules  where  available 
for  sampling). 

A  significant  conclusion  from  these  results  is  the  smallness  of  the  impact  of  flow  transience  and  stage 
motion  on  the  gas  flow  in  the  plume-atmosphere  mixing  layer.  Aside  from  statistical  scatter  that  is  visible 
in  the  three  transient  solutions,  there  is  practically  no  difference  in  the  location  and  properties  of  the  mixing 
layer  among  the  four  cases.  Note  the  moving  stage  result  has  larger  statistical  scatter  than  the  plume  only 
and  fixed  body  cases,  since  the  macroparameters  were  sampled  over  0.1  ms  time  in  the  former  case  and 


0.5  ms  in  the  latter  ones.  The  time  to  reach  steady  state  for  the  first  50  m  of  the  mixing  layer  is  on  the  order 


of  0.01  s,  which 


explains  the  similarity  of  the  results.  A  smaller  free  stream  velocity  would  increase  the  time 


to  reach  steady 


state,  but  still  this  would  be  insignificant  compared  to  the  time  scale  of  0.3  s. 
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Figure  3.  Liquid  propellant  thruster:  plume  H2O  number  density  fields  (mol/m3)  at  0.3  s  after  thruster 
ignition. 


Consider  now  the  impact  of  the  stage  motion  on  the  plume  species  density  distributions.  The  number 
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density  fields  of  H2O  at  0.3  s  after  plume  ignition  are  shown  in  Fig.  3.  There  is  some  interaction  of  the 
plume-atmosplrere  mixing  layer  and  the  weak  shock  generated  by  the  plume-stage  impingement  interaction 
that  results  in  a  small  increase  in  the  number  density  at  about  50  m  from  the  nozzle.  This  increase  is  not 
clearly  visible  in  the  moving  stage  case,  where  the  DSMC  statistical  noise  is  too  high  for  subtle  details  of 
the  flow  to  be  discernible.  Generally,  the  structure  of  the  plume-atmosphere  and  plume-stage  impingement 
shocks  and  the  near  field  of  the  plume  is  quantitatively  similar  in  the  transient  and  steady  state  cases.  There 


are  some  vortex-like  structures  in  the  mixing  layer  noticeable  in  both  carbon  monoxide  and  water  transient 


solutions,  but  the  nature  and  details  of  these  structures  are  out  of  the  scope  of  this  work. 
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Figure  4.  Liquid  propellant  thruster:  translational  temperature  fields  (K)  at  1.0  s  after  thruster  ignition. 

Due  to  accelerations  on  the  lower  stage,  the  lower  stage  relative  velocity  increases  significantly  with  time. 
Starting  from  2  m/s,  the  velocity  increases  to  about  25  m/s  at  0.3  s  and  then  to  about  75  m/s  at  1  s.  This 
relatively  high  velocity  is  still  not  sufficient  to  affect  the  gas  properties  during  stage  separation.  This  is 
illustrated  in  Fig.  4  where  the  translational  temperature  fields  are  given  for  the  moving  stage  and  steady 
state  cases  at  1  second.  The  other  two  cases  considered  above  are  not  shown  here  since  they  are  similar  to 
the  ones  shown.  The  structure  of  the  freestream-plume  interaction  region  is  very  similar  for  the  transient 
and  steady-state  cases,  and  resembles  strongly  that  of  the  early  time  moment  of  0.33  s.  The  results  indicate 
that  the  steady  state  result  captures  all  main  features  of  the  transient  one,  whereas  the  transient  case  fails 
to  adequately  represent  the  flow  near  the  lower  stage  due  to  poor  statistics.  The  total  number  of  particles 
of  5  million  used  in  the  computations  and  the  average  number  of  particles  per  sampling  cell  of  about  one 
is  obviously  too  small  for  a  quality  temperature  evaluation.  For  the  transient  case,  it  is  clear  that  a  three- 

11  of  21 


American  Institute  of  Aeronautics  and  Astronautics 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


dimensional  modeling  would  greatly  amplify  this  problem,  making  accurate  temperature  estimation  near 


the  lower  stage  at  higher  separation  distances  prohibitively  expensive  from  the  computational  standpoint  at 


longer  times. 
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Figure  5.  Liquid  propellant  thruster:  vibrational  temperature  field  (K)  at  1.0  s  after  thruster  ignition. 


The  statistical  fluctuations  are  even  more  severe  when  vibrational  temperature  is  considered.  Accurate 
prediction  of  vibrational  temperatures  is  crucial  for  analysis  of  infrared  (IR)  signatures.  The  computations 
conducted  for  a  discrete  model  of  the  molecular  internal  degrees  of  freedom  where  we  include  many  levels 
(not  shown  here)  indicate  that  the  vibrational  temperature  fluctuations  may  be  too  strong  for  any  meaningful 
analysis  to  be  conducted.  For  the  continuous  internal  energy  model,  the  actual  impact  of  statistical  scatter  is 
smaller  but  still  very  strong,  as  illustrated  in  Fig.  5  for  the  water  vibrational  temperature.  The  temperature 
was  averaged  over  0.05  ms  in  the  transient  case  in  order  to  reduce  the  statistical  noise.  Still,  the  mixing 
layer  and  the  lower  stage  region  temperatures  are  noticeably  smaller  than  in  the  steady  state  solution.  The 
impact  of  transience  on  radiation  signatures  may  not  be  that  significant  since  the  maximum  difference  is 
observed  in  regions  with  very  low  densities  of  H20. 

As  a  step  toward  a  complete  treatment  of  the  radiation  environment  we  have  performed  calculations 
with  the  Air  Force  chemistry,  radiation  and  signature  code  SOCRATES-P.11  In  Figure  6  we  show  broadside 
infrared  images  of  radiance  from  vibrationally  excited  gases  as  could  be  viewed  from  a  far  away  sensor 
computed  with  SOCRATES-P.  The  calculation  was  performed  with  the  same  parameters  as  the  steady-state 
SMILE  calculation,  but  using  a  smaller  simulation  volume  in  quarter-symmetry,  using  1.2  million  cells  in  a 
regionally-adapted  grid  populated  by  13  million  molecules.  The  vibrational  excitations  are  described  with 
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c 

Figure  6.  Infrared  images  of  the  liquid-propellant  staging,  as  computed  with  SOCRATES-P  for  4.7-pm  excited 
CO  (left)  and  6.3-pm  excited  I  i  O  (right).  The  width  of  each  image  is  30  m.  The  intensity  scale  has  been  dis¬ 
torted  from  linear  using  an  image-processing  gamma=0.01  transfer-function  to  show  the  radiance  enhancement 
near  the  lower  stage. 

two-state  (ground  and  vibrationally  excited)  models,  where  collisional  excitation,  radiation  and  quenching 
are  included.  In  either  band,  the  dominant  intensity  source  is  the  core  flow  exiting  the  upper-stage  engine,  but 
the  H2O  radiance  shows  the  shock  structure  more  clearly.  The  slight  chevron-pattern  in  the  water  radiance 
imagery  is  an  artifact  of  the  plume-atmosphere  shock  leaving  the  solution  domain,  as  viewed  broadside. 
These  SOCRATES-P  calculations  are  the  first  step  in  examining  the  radiation  environment  and  signature 
for  staging  scenarios. 


Figure  7.  Liquid  propellant  thruster:  temporal  profiles  of  stage  separation  distance  and  plume  force  on  lower 
stage. 


The  computations  have  shown  that  the  lower  stage  trajectory  is  mostly  governed  by  the  thrust  force  T, 
with  a  smaller  contribution  of  the  plume  force  Fp.  The  distance  between  the  upper  and  lower  stages  as 
a  function  of  time  is  plotted  in  Fig.  7  for  the  actual  stage  separation  observed  in  the  computations  and  a 
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separation  calculated  using  only  the  initial  stage  velocity  and  thrust  force  propelling  the  upper  stage  (no 
plume  force).  The  ignition  starts  at  0.02s,  and  at  1  s  the  distance  between  the  stages  is  by  about  15%  larger 
when  the  plume  force  is  included.  The  plume  force  has  the  biggest  impact  on  the  stage  velocity  during  the 
first  0.1  s,  when  the  separation  distance  is  still  small  and  the  plume  force  impingement  to  lower  stage  mass 
ratio  is  comparable  to  the  thrust  to  upper  stage  mass  ratio.  Note  that  even  immediately  after  firing,  the 
plume  force  on  the  lower  stage  is  only  about  20%  of  the  total  thrust  force,  both  due  to  the  plume  divergence 
after  the  nozzle  exit  and  force  decrease  in  the  shock  front. 

An  important  conclusion  from  the  modeling  results  of  a  25  kN  liquid  propellant  thruster  stage  separation 
is  that  the  impact  of  the  stage  motion  and  unsteady  effects  on  flow  properties,  both  close  and  far  from 
the  moving  stage,  is  small.  This  indicates  that  quasi-steady  state  treatments  for  such  scenarios  may  be 
quite  adequate.  Moreover,  potential  benefits  of  accurate  prediction  of  the  stage  trajectory  with  transient 
DSMC  modeling  may  be  hindered  by  difficulties  of  flowheld  and  signature  predictions  associated  with  DSMC 
statistical  scatter  in  transient  flow  modeling  with  the  DSMC  method.  We  also  find  that  the  impact  of  the 
plume  force  on  the  stage  separation  is  noticeable  and  generally  can  not  be  neglected. 

VI.  Impact  of  Statistical  Scatter  on  Stage  Trajectory 

The  application  of  the  DSMC  method  to  modeling  stage  separation  may  imply  a  significant  level  of 
statistical  scatter  in  the  computation  of  the  force  from  the  plume  impingement  on  the  lower  stage.  In  the 
computations  presented  in  the  previous  section,  the  plume  force  was  estimated  based  on  no  more  than  a 
few  thousand  gas-surface  collisions  even  for  smaller  separation  distances.  For  larger  distances,  the  number 
of  surface  collisions  decreases  by  almost  three  orders  of  magnitude,  although  the  impact  of  the  plume  force 
compared  to  the  thrust  force  also  decreases  drastically.  Generally,  a  strong  dependence  of  the  stage  trajectory 
on  the  accuracy  of  the  force  calculation  signals  that  the  number  of  particles  used  is  not  adequate,  and  needs 
to  be  increase,  or  a  spatial  weighting  approach  to  increase  this  accuracy  may  have  to  be  developed  and 
applied. 

To  analyze  the  possible  contribution  of  statistical  scatter  to  the  stage  trajectory  predictions,  the  DSMC 
computations  have  been  performed  for  different  numbers  of  simulated  particles.  This  set  of  computations 
have  been  performed  using  the  same  nozzle  exit  flow  properties  as  for  the  above  axisymmetric  computations, 
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but  in  a  two-dimensional  Cartesian  treatment.  A  2D  Cartesian  consideration  allows  one  to  include  body 


rotation  as  well  as  deviation  from  the  plume  axis,  not  possible  in  an  axisymmetric  modeling.  Note  that  the 
use  of  the  same  number  density  and  flow  velocity  at  the  nozzle  exit  results  in  an  effectively  much  higher 
thrust  at  the  lower  stage  in  2D,  thus  amplifying  the  effect  of  the  plume. 
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Figure  8.  2D  liquid  propellant  case:  the  impact  of  the  number  of  simulated  particles  on  stage  separation 
distance  (left),  deviation  from  the  plume  axis  (center),  and  the  rotation  angle  (right). 

The  computations  are  conducted  for  three  values  of  the  total  number  of  simulated  molecules,  5  million, 
1  million,  and  0.2  million.  Analysis  of  the  results  have  shown  that  gas  flowfield  solutions  are  qualitatively 
similar  for  these  three  cases,  while  a  quantitative  comparison  is  complicated  by  significant  statistical  scatter 
when  lower  numbers  of  particles  are  used.  The  influence  of  the  scatter  in  the  plume  force  evaluation  on  the 
total  separation  distance  measured  as  the  closest  distance  between  the  two  stages  is  given  in  Fig.  8(left). 
Although  there  is  a  visible  difference  between  the  separation  distance  for  different  numbers  of  particles,  these 
differences  are  much  smaller  than  their  difference  from  the  case  when  the  stage  moved  only  due  to  the  thrust 
propelling  the  upper  stage. 

In  an  ideal  case  of  a  fully  symmetric  plume,  there  should  be  no  deviation  of  the  stage  center  of  mass  from 
the  plume  axis.  The  moving  lower  stage  should  always  be  aligned  with  this  axis  since  the  lift  force,  as  well 
as  the  pitching  moment,  are  zero.  However,  the  statistical  noise  inherent  in  any  DSMC  modeling  result  in 
the  lift  and  torque  forces  on  the  moving  stage.  Moreover,  the  stage  translation  off  the  axis  as  well  as  the 
rotation  further  increase  these  forces  due  to  created  asymmetry.  The  effect  of  the  statistical  scatter  on  the 
distance  between  the  lower  stage  center  of  mass  and  the  plume  axis  is  given  in  Fig.  8(center).  The  conclusion 
here  is  that  the  impact  of  the  number  of  particles  on  the  stage  deviation  from  the  axis  is  negligible.  This 
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is  because  the  lift  force  is  small  compared  to  the  drag  force,  and  the  positive  and  negative  values  of  the  lift 
force  at  consecutive  displacements  of  the  lower  stage  nearly  equilibrate  over  time. 

Similarly,  the  rotation  of  the  stage  is  insignificant  over  the  first  second  after  stage  ignition,  as  illustrated 
in  Fig.  8(right).  It  is  also  important  to  mention  that  the  results  of  this  section  are  2D,  and  the  impact  of 
the  plume  force  is  strongly  amplified,  by  a  factor  of  five,  as  compared  to  the  corresponding  axisymmetric  - 
and  3D  -  case.  Such  a  weak  dependence  of  the  lower  stage  trajectory  on  statistical  fluctuations  caused  by 
limited  number  of  gas-surface  collisions  is  an  indication  that  even  a  full  3D  DSMC  modeling  is  possible  for 
a  moving  lower  stage  with  an  acceptable  accuracy  of  trajectory  estimation. 

VII.  Solid  propellant  thrusters 

The  weak  dependence  of  the  stage  separation  modeling  results  on  the  stage  motion  and  unsteady  nature 
of  the  flow  is  primarily  related  to  the  relatively  fast  time  to  reach  steady  state  both  in  the  plume-stage  and 
the  plume-free  stream  interaction  regions.  Usually,  the  flow  needs  a  longer  time  to  reach  steady  state  for 
solid  propellant  rocket  thrusters,  since  the  plume  is  two-phase  and  includes  both  gas  and  alumina  particle 
transport.  The  particulates  typically  move  slower  than  gas,  and  there  is  also  some  gas-particle  interaction, 
which  both  result  in  longer  transient  time  than  for  liquid  propellant  thrusters  under  similar  conditions. 
To  study  the  effect  of  particulates  on  the  transient  stage  separation,  the  DSMC  computations  have  been 
conducted  for  a  35  kN  solid  propellant  thruster  with  stagnation  gas  parameters  similar  to  the  previously 
considered  liquid  propellant  case,  and  two  sizes  of  particulates,  small  and  large  with  diameters  of  0.2  /im 
and  3.6  /xm  diameter,  respectively.  As  mentioned  above,  a  close  to  maximum  particle  loading  of  20%  was 
used  in  order  to  emphasize  the  effect  of  particles. 

The  gas  and  particle  properties  at  the  nozzle  exit  plane  obtained  with  CFD+- 1-  continuum  solver  are 
shown  in  Fig.  9.  As  gas  density  profile  clearly  shows,  there  is  a  noticeable  boundary  layer  formed  at  the 
nozzle  surface.  The  density  of  the  small  particulates  is  qualitatively  very  similar  to  the  gas  density,  which 
is  because  the  mass  of  these  particles  is  too  small  to  be  different  from  the  gas  (the  gas  and  0.2  /xm  particle 
velocities  practically  coincide).  The  drag  force  on  larger  alumina  particles  is  not  high  enough  for  them  to 
move  along  gas  streamlines,  and  their  number  density  shape  differs  from  that  of  the  gas.  Note  that  only  the 
large  particle  surface  temperature  is  presented  in  Fig.  9,  since  the  surface  temperature  of  smaller  particles 
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2000 


Figure  9.  Solid  propellant  thruster:  flow  properties  across  the  nozzle  exit  plane. 


coincide  with  that  of  the  gas.  It  is  important  to  note  that  both  the  small  and  large  particle  temperature 
is  noticeably  lower  than  the  alumina  melting  temperature  of  2,325  K,  which  means  that  the  particles  have 
already  solidified  traveling  in  the  diverging  part  of  the  nozzle. 
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Figure  10.  Solid  propellant  thruster:  gas  translational  temperature  (K)  at  0.3  s  after  thruster  ignition.  The 


top  half  are  the  unsteady  results  and  the  bottom  half  are  the  steady-state  results. 


We  now  consider  the  effect  of  the  stage  motion  and  unsteady  effects  on  gas  and  particulate  properties. 
The  gas  temperature  fields  for  a  stage  separation  at  0.3  s  after  the  thruster  ignition  are  shown  in  Fig.  10. 
The  distance  between  the  stage  at  this  time  moment  is  about  7  m.  Here  and  below,  the  top  half  illustrates 
the  results  of  the  transient  flow  development  with  the  moving  lower  stage,  and  the  lower  half  shown  the 
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results  of  the  steady  state  modeling  with  the  distance  between  the  stages  that  corresponds  to  the  transient 


case.  Comparison  of  the  results  shows  that  there  is  no  significant  impact  of  the  flow  transience,  similar  to 


the  liquid  case  considered  above.  The  gas  temperature  in  the  wake  is  somewhat  lower  in  the  transient  case, 


which  is  attributed  to  the  ensemble-averaging  sampling  over  a  relatively  small  number  of  particles  in  this 


rarefied  region  and  not  the  actual  motion  of  the  stage. 
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Figure  11.  Solid  propellant:  number  density  fields  (molecule/m3)  of  large  particulates  at  0.3  seconds  after 
separation.  The  top  half  is  the  unsteady  results  and  the  lower  half  is  the  steady  state  results. 


Comparison  of  the  transient  and  steady  state  number  density  fields  of  large  particles  at  0.3  s  is  given  in 
Fig.  11.  It  is  clear  that  generally  there  is  a  very  good  agreement  between  the  two  solutions.  Note  that  the 
results  for  smaller  particles  agree  as  well.  For  both  particle  sizes  there  is  an  increase  in  the  number  density  in 
the  shock  front  near  the  lower  stage.  Smaller  particles  follow  the  gas  streamlines  in  the  bow  shock  and  then 
in  the  expansion  region.  Due  to  their  large  mass,  the  3.6/xrn  particles  reflected  from  the  body  are  impacted 
by  the  gas  to  a  lesser  extent  than  the  smaller  ones,  and  do  not  follow  gas  streamlines  in  the  bow  shock.  The 
trajectories  of  the  large  particles  that  did  not  collide  with  the  stage  are  weakly  changed  by  the  bow  shock. 
Note  also  that  the  particle  temperatures  do  not  change  significantly  throughout  the  computational  domain 
outside  the  shock  region  in  front  of  the  lower  stage.  In  that  small  region,  the  smaller  particle  temperature 
first  increases  above  the  melting  temperatures,  and  then  drops  to  below  melting  temperatures  near  the 
surface.  This  indicates  that  the  smaller  particles  are  likely  to  melt  and  be  in  the  liquid  states  during  their 
collisions  with  the  body.  The  particle  melting  was  not  modeled  in  this  work.  The  temperature  of  the  larger 
particles  reaches  a  maximum  of  about  2280  K  inside  the  shock.  Generally,  special  care  needs  to  be  taken 
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when  modeling  the  interaction  of  particles  traveling  through  the  shock  with  the  stage  surface,  but  the  details 
of  this  process  were  not  in  the  scope  of  this  work.  The  important  conclusion  here  is  that  the  gas  and  particle 
behavior  is  not  visibly  impacted  by  the  stage  motion. 

VIII.  Conclusions 

The  direct  simulation  Monte  Carlo  method  has  been  used  to  simulate  the  staging  process  of  a  generic 
sounding  rocket  at  an  altitude  of  100  km.  Both  liquid  and  solid  propellants  have  been  modeled,  with  thrust 
levels  of  25  kN  and  34  kN,  respectively.  A  CFD++  commercial  Navier-Stokes  solver  was  used  to  predict  the 
gas  and  particulate  properties  at  the  nozzle  exit  plane.  The  2D/axisynnnetric  module  of  the  DSMC  based 
solver  SMILE  was  then  used  to  calculate  the  flow  from  the  exit  plane.  For  these  calculations,  SMILE  was 
extended  to  include  the  body  motion,  and  all  unsteady  flow  effects  were  included.  Both  constant  thrust  force 
and  time-dependent  plume  force  on  the  lower  stage  were  considered.  The  acceleration  of  the  lower  stage 
relative  to  the  upper  stage  varied  from  about  \Ag  immediately  after  the  ignition  to  about  7g  at  1  s  after 
ignition. 

From  our  DSMC  stage  separation  results,  we  have  seen  that  the  engine  start-up  produces  a  plume  that 
establishes  a  shock  with  the  atmosphere.  The  plume’s  impingement  on  the  lower  stage  creates  a  second 
shock  and  a  stagnation  region  near  the  plume  impact  on  the  lower  stage.  This  second  shock  is  established 
quickly,  and  it  persists  but  becomes  weaker  as  the  stages  move  away  from  each  other.  Although  the  two 
shock  structures  may  intersect  weakly  downstream  (see  Fig.  3  at  40  m  downstream  of  the  nozzle  exit),  there 
is  almost  no  interaction  between  the  two  shock  structures.  To  analyze  the  impact  of  the  lower  stage  motion, 
in  addition  to  the  fully  transient  case,  computations  have  been  performed  for  a  transient  plume  with  a  fixed 
body  within  a  steady  state  treatment  (i.e.  quasi-steady  assumptions),  with  the  stage  separation  distance 
corresponding  to  that  of  the  fully  transient  case.  The  computations  showed  that  the  stage  motion  and  other 
unsteady  flow  effects  have  little  influence  on  the  results  of  the  computations  as  compared  to  the  steady  state 
cases.  This  result  is  consistent  with  the  much  faster  time-scale  of  the  molecular  gas-gas  collisions  than  the 
time-scale  for  motion  of  the  solid  body.  This  indicates  that  quasi-steady  state  modeling  may  be  preferential 
due  to  the  small  impact  of  statistical  scatter  for  these  scenarios.  There  are  some  small  differences  between 
the  steady  and  unsteady  solutions  in  flow  regions  with  relatively  small  numbers  of  simulated  particles.  These 
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differences  are  attributed  to  the  poorer  statistics  of  the  transient  solutions  as  compared  to  the  steady  state 
solutions.  These  conclusions  hold  both  for  liquid  and  solid  propellant  rockets,  and  may  be  generalized  to 
staging  processes  with  lower  stage  accelerations. 

It  is  important  to  establish  the  need  to  explicitly  include  time-dependence  in  such  fluid  flow  calculations, 
as  it  can  incur  a  large  computational  expense.  One  approximate  way  to  do  this  before  undertaking  potentially 
time-consuming  computations  is  to  examine  the  basic  time-scales  of  motion.  A  fundamental  time  scale  for 
gas-gas  collisions  is  the  local  molecular  collision  time  t.  If  the  solid  body  can  move  a  substantial  distance 
within  a  molecular  collision  period,  the  motion  of  the  body  will  be  coupled  to  the  gas  flow.  A  way  to  quantify 
the  relative  time  scales  of  the  molecular  collision  period  and  the  solid  body  motion  is  the  ratio  of  the  mean 
free  path,  A,  divided  by  the  collision  time,  t,  to  the  speed  of  the  solid  body,  vSoiid-  A /{t  ■  vsoud)-  When  this 
ratio  is  much  greater  than  one  the  gas-gas  collision  times  are  much  less  than  the  solid  body  motion,  and  in 
this  case  a  quasi-steady  state  treatment,  where  the  solid  body  motion  is  “frozen”  along  its  trajectory  and 
several  independent  steady  state  calculations  are  done,  should  be  sufficient.  For  the  present  scenarios  with 
the  number  density  on  the  order  of  5  x  1020  m  ,  a  mean  molecular  speed  ~1500  m/s,  and  using  a  molecular 
cross  section  of  10”19  nr2,  we  find  a  ratio  of  ~  60  ^  1.  This  confirms  that  a  quasi-steady  state  treatment 
should  be  nearly  equivalent  to  a  fully  coupled  time-dependent  treatment  of  the  body  motion,  as  we  observe. 
Of  course,  this  simplified  approach  is  approximate  since  it  does  not  take  into  consideration  the  actual  scales 
of  the  flow  development  that  generally  should  be  accounted  for. 

The  impact  of  the  statistical  fluctuations  on  the  stage  separation  trajectory  has  been  studied  for  a  2D  flow, 
and  was  found  relatively  small  compared  to  the  general  impact  of  the  plume  force  on  the  stage  separation 
distance.  This  conclusion  is  applicable  to  separation  distance,  motion  in  the  transverse  direction,  and 
rotation.  These  findings  are  explained  by  the  fact  that  the  motion  due  to  the  lift  force  and  pitching  moment 
is  small  compared  to  that  caused  by  the  drag  force.  The  positive  and  negative  values  of  the  lift  force  and 
pitching  moment  at  consecutive  displacements  of  the  lower  stage  due  to  statistical  scatter  nearly  equilibrate 
over  time.  These  results  may  also  be  an  indication  that  similar  transient  scenarios  can  be  examined  in  full 
3D  with  little  concern  over  DSMC  statistical  errors  from  surface  force  fluctuations. 

As  a  step  toward  a  complete  treatment  of  the  staging  radiation  environment,  we  have  also  performed 
calculations  with  the  Air  Force  chemistry,  radiation  and  signature  code  SOCRATES-P.11  IR  radiance  imagery 
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generated  from  these  preliminary  calculations  reveal  the  plume-atmosphere  shock  structure  and  plume-lower 
stage  impingement  radiation.  These  SOCRATES-P  calculations  are  the  first  step  in  examining  the  full 
transient  radiation  environment  and  signature  for  staging  scenarios. 
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